{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# surface_energy_static calculation style\n",
    "\n",
    "**Lucas M. Hale**, [lucas.hale@nist.gov](mailto:lucas.hale@nist.gov?Subject=ipr-demo), *Materials Science and Engineering Division, NIST*.\n",
    "\n",
    "## Introduction\n",
    "\n",
    "The surface_energy_static calculation style evaluates the formation energy for a free surface by slicing an atomic system along a specific plane.\n",
    "\n",
    "### Version notes\n",
    "\n",
    "- 2018-07-09: Notebook added.\n",
    "- 2019-07-30: Description updated and small changes due to iprPy version.\n",
    "- 2020-05-22: Version 0.10 update - potentials now loaded from database.\n",
    "- 2020-09-22: Calculation updated to use atomman.defect.FreeSurface class. Setup and parameter definition streamlined.\n",
    "\n",
    "### Additional dependencies\n",
    "\n",
    "### Disclaimers\n",
    "\n",
    "- [NIST disclaimers](http://www.nist.gov/public_affairs/disclaimer.cfm)\n",
    "- Other atomic configurations at the free surface for certain planar cuts may have lower energies.  The atomic relaxation will find a local minimum, which may not be the global minimum.  Additionally, the material cut is planar perfect and therefore does not explore the effects of atomic roughness.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Method and Theory\n",
    "\n",
    "First, an initial system is generated.  This is accomplished by\n",
    "\n",
    "1. Starting with a unit cell system.\n",
    "\n",
    "2. Generating a transformed system by rotating the unit cell such that the new system's box vectors correspond to crystallographic directions, and filled in with atoms to remain a perfect bulk cell when the three boundaries are periodic.\n",
    "\n",
    "3. All atoms are shifted by a fractional amount of the box vectors if needed.\n",
    "\n",
    "4. A supercell system is constructed by combining multiple replicas of the transformed system.\n",
    "\n",
    "Two LAMMPS simulations are then performed that apply an energy/force minimization on the system, and the total energy of the system after relaxing is measured, $E_{total}$.  In the first simulation, all of the box's directions are kept periodic (ppp), while in the second simulation two are periodic and one is non-periodic (ppm).  This effectively slices the system along the boundary plane creating two free surfaces, each with surface area\n",
    "\n",
    "$$A = \\left| \\vec{a_1} \\times \\vec{a_2} \\right|,$$\n",
    "\n",
    "where $\\vec{a_1}$ and $\\vec{a_2}$ are the two lattice vectors corresponding to the periodic in-plane directions.\n",
    "\n",
    "The formation energy of the free surface, $E_{f}^{surf}$, is computed in units of energy over area as\n",
    "\n",
    "$$E_{f}^{surf} = \\frac{E_{total}^{surf} - E_{total}^{0}} {2 A}.$$\n",
    "\n",
    "The calculation method allows for the specification of which of the three box dimensions the cut is made along.  If not specified, the default behavior is to make the $\\vec{c}$ vector direction non-periodic.  This choice is due to the limitations of how LAMMPS defines triclinic boxes.  $\\vec{c}$ vector is the only box vector that is allowed to have a component in the Cartesian z direction.  Because of this, the other two box vectors are normal to the z-axis and therefore will be in the cut plane.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Demonstration"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1. Setup"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.1. Library imports"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Import libraries needed by the calculation. The external libraries used are:\n",
    "\n",
    "- [numpy](http://www.numpy.org/)\n",
    "\n",
    "- [atomman](https://github.com/usnistgov/atomman)\n",
    "\n",
    "- [iprPy](https://github.com/usnistgov/iprPy)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Notebook last executed on 2020-09-22 using iprPy version 0.10.2\n"
     ]
    }
   ],
   "source": [
    "# Standard library imports\n",
    "from pathlib import Path\n",
    "import os\n",
    "import shutil\n",
    "import datetime\n",
    "\n",
    "# http://www.numpy.org/\n",
    "import numpy as np \n",
    "\n",
    "# https://github.com/usnistgov/atomman \n",
    "import atomman as am\n",
    "import atomman.lammps as lmp\n",
    "import atomman.unitconvert as uc\n",
    "\n",
    "# https://github.com/usnistgov/iprPy\n",
    "import iprPy\n",
    "\n",
    "print('Notebook last executed on', datetime.date.today(), 'using iprPy version', iprPy.__version__)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 1.2. Default calculation setup"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Specify calculation style\n",
    "calc_style = 'surface_energy_static'\n",
    "\n",
    "# If workingdir is already set, then do nothing (already in correct folder)\n",
    "try:\n",
    "    workingdir = workingdir\n",
    "\n",
    "# Change to workingdir if not already there\n",
    "except:\n",
    "    workingdir = Path('calculationfiles', calc_style)\n",
    "    if not workingdir.is_dir():\n",
    "        workingdir.mkdir(parents=True)\n",
    "    os.chdir(workingdir)\n",
    "    \n",
    "# Initialize connection to library\n",
    "with open('C:/Users/lmh1/Documents/potentials_nist_gov/password.txt') as f:\n",
    "    user, pswd = f.read().strip().split()\n",
    "library = iprPy.Library(load=['lammps_potentials'], username=user, password=pswd)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2. Assign values for the calculation's run parameters"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 2.1. Specify system-specific paths\n",
    "\n",
    "- __lammps_command__ is the LAMMPS command to use (required).\n",
    "\n",
    "- __mpi_command__ MPI command for running LAMMPS in parallel. A value of None will run simulations serially."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "lammps_command = 'lmp_serial'\n",
    "mpi_command = None"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 2.2. Load interatomic potential\n",
    "\n",
    "- __potential_name__ gives the name of the potential_LAMMPS reference record in the iprPy library to use for the calculation.  \n",
    "\n",
    "- __potential__ is an atomman.lammps.Potential object (required)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "potential_name = '1999--Mishin-Y--Ni--LAMMPS--ipr1'\n",
    "\n",
    "# Retrieve potential and parameter file(s)\n",
    "potential = library.get_lammps_potential(id=potential_name, getfiles=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 2.3. Load initial unit cell system\n",
    "\n",
    "- __ucell__ is an atomman.System representing a fundamental unit cell of the system (required).  Here, this is loaded from the database for the prototype."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "avect =  [ 3.520,  0.000,  0.000]\n",
      "bvect =  [ 0.000,  3.520,  0.000]\n",
      "cvect =  [ 0.000,  0.000,  3.520]\n",
      "origin = [ 0.000,  0.000,  0.000]\n",
      "natoms = 4\n",
      "natypes = 1\n",
      "symbols = ('Ni',)\n",
      "pbc = [ True  True  True]\n",
      "per-atom properties = ['atype', 'pos']\n",
      "     id |   atype |  pos[0] |  pos[1] |  pos[2]\n",
      "      0 |       1 |   0.000 |   0.000 |   0.000\n",
      "      1 |       1 |   0.000 |   1.760 |   1.760\n",
      "      2 |       1 |   1.760 |   0.000 |   1.760\n",
      "      3 |       1 |   1.760 |   1.760 |   0.000\n"
     ]
    }
   ],
   "source": [
    "# Create ucell by loading prototype record\n",
    "ucell = am.load('crystal', potential=potential, family='A1--Cu--fcc', database=library)\n",
    "\n",
    "print(ucell)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 2.4. Specify the defect parameters\n",
    "\n",
    "- __hkl__ gives the Miller (hkl) or Miller-Bravais (hkil) plane to create the free surface on.\n",
    "\n",
    "- __cutboxvector__ specifies which of the three box vectors ('a', 'b', or 'c') is to be made non-periodic to create the free surface. \n",
    "\n",
    "- __shiftindex__ can be used for complex crystals to specify different termination planes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "hkl = [1,0,0]\n",
    "cutboxvector = 'c'\n",
    "shiftindex = 0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 2.5. Modify system\n",
    "\n",
    "- __sizemults__  list of three integers specifying how many times the ucell vectors of $a$, $b$ and $c$ are replicated in creating system.\n",
    "- __minwidth__ specifies a minimum width that the system should be along the cutboxvector direction. The given sizemult in that direction will be increased if needed to ensure that the system is at least this wide.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "sizemults = [5, 5, 10]\n",
    "minwidth = uc.set_in_units(0.0, 'angstrom')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 2.6. Specify calculation-specific run parameters\n",
    "\n",
    "- __energytolerance__ is the energy tolerance to use during the minimizations. This is unitless.\n",
    "\n",
    "- __forcetolerance__ is the force tolerance to use during the minimizations. This is in energy/length units.\n",
    "\n",
    "- __maxiterations__ is the maximum number of minimization iterations to use.\n",
    "\n",
    "- __maxevaluations__ is the maximum number of minimization evaluations to use.\n",
    "\n",
    "- __maxatommotion__ is the largest distance that an atom is allowed to move during a minimization iteration. This is in length units."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "energytolerance = 1e-8\n",
    "forcetolerance = uc.set_in_units(0.0, 'eV/angstrom')\n",
    "maxiterations = 10000\n",
    "maxevaluations = 100000\n",
    "maxatommotion = uc.set_in_units(0.01, 'angstrom')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3. Define calculation function(s) and generate template LAMMPS script(s)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 3.1. min.template"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "with open('min.template', 'w') as f:\n",
    "    f.write(\"\"\"#LAMMPS input script that performs an energy minimization\n",
    "\n",
    "box tilt large\n",
    "\n",
    "<atomman_system_pair_info>\n",
    "\n",
    "thermo_style custom step lx ly lz pxx pyy pzz pe\n",
    "thermo_modify format float %.13e\n",
    "\n",
    "compute peatom all pe/atom \n",
    "\n",
    "min_modify dmax <dmax>\n",
    "\n",
    "dump dumpit all custom <maxeval> atom.* id type x y z c_peatom\n",
    "dump_modify dumpit format <dump_modify_format>\n",
    "\n",
    "minimize <etol> <ftol> <maxiter> <maxeval>\"\"\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 3.2. surface_energy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "def surface_energy(lammps_command, ucell, potential, hkl,\n",
    "                   mpi_command=None, sizemults=None, minwidth=None, even=False,\n",
    "                   conventional_setting='p', cutboxvector='c', \n",
    "                   atomshift=None, shiftindex=None,\n",
    "                   etol=0.0, ftol=0.0, maxiter=10000,\n",
    "                   maxeval=100000, dmax=uc.set_in_units(0.01, 'angstrom')):\n",
    "    \"\"\"\n",
    "    Evaluates surface formation energies by slicing along one periodic\n",
    "    boundary of a bulk system.\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    lammps_command :str\n",
    "        Command for running LAMMPS.\n",
    "    ucell : atomman.System\n",
    "        The crystal unit cell to use as the basis of the stacking fault\n",
    "        configurations.\n",
    "    potential : atomman.lammps.Potential\n",
    "        The LAMMPS implemented potential to use.\n",
    "    hkl : array-like object\n",
    "        The Miller(-Bravais) crystal fault plane relative to ucell.\n",
    "    mpi_command : str, optional\n",
    "        The MPI command for running LAMMPS in parallel.  If not given, LAMMPS\n",
    "        will run serially.\n",
    "    sizemults : list or tuple, optional\n",
    "        The three System.supersize multipliers [a_mult, b_mult, c_mult] to use on the\n",
    "        rotated cell to build the final system. Note that the cutboxvector sizemult\n",
    "        must be an integer and not a tuple.  Default value is [1, 1, 1].\n",
    "    minwidth : float, optional\n",
    "        If given, the sizemult along the cutboxvector will be selected such that the\n",
    "        width of the resulting final system in that direction will be at least this\n",
    "        value. If both sizemults and minwidth are given, then the larger of the two\n",
    "        in the cutboxvector direction will be used. \n",
    "    even : bool, optional\n",
    "        A True value means that the sizemult for cutboxvector will be made an even\n",
    "        number by adding 1 if it is odd.  Default value is False.\n",
    "    conventional_setting : str, optional\n",
    "        Allows for rotations of a primitive unit cell to be determined from\n",
    "        (hkl) indices specified relative to a conventional unit cell.  Allowed\n",
    "        settings: 'p' for primitive (no conversion), 'f' for face-centered,\n",
    "        'i' for body-centered, and 'a', 'b', or 'c' for side-centered.  Default\n",
    "        behavior is to perform no conversion, i.e. take (hkl) relative to the\n",
    "        given ucell.\n",
    "    cutboxvector : str, optional\n",
    "        Indicates which of the three system box vectors, 'a', 'b', or 'c', to\n",
    "        cut with a non-periodic boundary (default is 'c').\n",
    "    atomshift : array-like object, optional\n",
    "        A Cartesian vector shift to apply to all atoms.  Can be used to shift\n",
    "        atoms perpendicular to the fault plane to allow different termination\n",
    "        planes to be cut.  Cannot be given with shiftindex.\n",
    "    shiftindex : int, optional\n",
    "        Allows for selection of different termination planes based on the\n",
    "        preferred shift values determined by the underlying fault generation.\n",
    "        Cannot be given with atomshift. If neither atomshift nor shiftindex\n",
    "        given, then shiftindex will be set to 0.\n",
    "    etol : float, optional\n",
    "        The energy tolerance for the structure minimization. This value is\n",
    "        unitless. (Default is 0.0).\n",
    "    ftol : float, optional\n",
    "        The force tolerance for the structure minimization. This value is in\n",
    "        units of force. (Default is 0.0).\n",
    "    maxiter : int, optional\n",
    "        The maximum number of minimization iterations to use (default is \n",
    "        10000).\n",
    "    maxeval : int, optional\n",
    "        The maximum number of minimization evaluations to use (default is \n",
    "        100000).\n",
    "    dmax : float, optional\n",
    "        The maximum distance in length units that any atom is allowed to relax\n",
    "        in any direction during a single minimization iteration (default is\n",
    "        0.01 Angstroms).\n",
    "    \n",
    "    Returns\n",
    "    -------\n",
    "    dict\n",
    "        Dictionary of results consisting of keys:\n",
    "        \n",
    "        - **'dumpfile_base'** (*str*) - The filename of the LAMMPS dump file\n",
    "          of the relaxed bulk system.\n",
    "        - **'dumpfile_surf'** (*str*) - The filename of the LAMMPS dump file\n",
    "          of the relaxed system containing the free surfaces.\n",
    "        - **'E_total_base'** (*float*) - The total potential energy of the\n",
    "          relaxed bulk system.\n",
    "        - **'E_total_surf'** (*float*) - The total potential energy of the\n",
    "          relaxed system containing the free surfaces.\n",
    "        - **'A_surf'** (*float*) - The area of the free surface.\n",
    "        - **'E_coh'** (*float*) - The cohesive energy of the relaxed bulk\n",
    "          system.\n",
    "        - **'E_surf_f'** (*float*) - The computed surface formation energy.\n",
    "    \n",
    "    Raises\n",
    "    ------\n",
    "    ValueError\n",
    "        For invalid cutboxvectors\n",
    "    \"\"\"\n",
    "    # Construct free surface configuration generator\n",
    "    surf_gen = am.defect.FreeSurface(hkl, ucell, cutboxvector=cutboxvector,\n",
    "                                     conventional_setting=conventional_setting)\n",
    "\n",
    "    # Check shift parameters\n",
    "    if shiftindex is not None:\n",
    "        assert atomshift is None, 'shiftindex and atomshift cannot both be given'\n",
    "        atomshift = surf_gen.shifts[shiftindex]\n",
    "    elif atomshift is None:\n",
    "        atomshift = surf_gen.shifts[0]\n",
    "\n",
    "    # Generate the free surface configuration\n",
    "    system = surf_gen.surface(shift=atomshift, minwidth=minwidth,\n",
    "                                  sizemults=sizemults, even=even)\n",
    "    A_surf= surf_gen.surfacearea\n",
    "\n",
    "    # Identify lammps_date version\n",
    "    lammps_date = lmp.checkversion(lammps_command)['date']\n",
    "\n",
    "    # Evaluate system with free surface\n",
    "    surf_results = relax_system(lammps_command, system, potential,\n",
    "                                mpi_command=mpi_command, etol=etol, ftol=ftol,\n",
    "                                maxiter=maxiter, maxeval=maxeval, dmax=dmax,\n",
    "                                lammps_date=lammps_date)\n",
    "    \n",
    "    # Extract results from system with free surface\n",
    "    dumpfile_surf = 'surface.dump'\n",
    "    shutil.move(surf_results['finaldumpfile'], dumpfile_surf)\n",
    "    shutil.move('log.lammps', 'surface-log.lammps')\n",
    "    E_total_surf = surf_results['potentialenergy']\n",
    "\n",
    "    # Evaluate perfect system (all pbc removes cut)\n",
    "    system.pbc = [True, True, True]\n",
    "    perf_results = relax_system(lammps_command, system, potential,\n",
    "                                mpi_command=mpi_command, etol=etol, ftol=ftol,\n",
    "                                maxiter=maxiter, maxeval=maxeval, dmax=dmax,\n",
    "                                lammps_date=lammps_date)\n",
    "    \n",
    "    # Extract results from perfect system\n",
    "    dumpfile_base = 'perfect.dump'\n",
    "    shutil.move(perf_results['finaldumpfile'], dumpfile_base)\n",
    "    shutil.move('log.lammps', 'perfect-log.lammps')\n",
    "    E_total_base = perf_results['potentialenergy']\n",
    "    \n",
    "    # Compute the free surface formation energy\n",
    "    E_surf_f = (E_total_surf - E_total_base) / (2 * A_surf)\n",
    "    \n",
    "    # Save values to results dictionary\n",
    "    results_dict = {}\n",
    "    \n",
    "    results_dict['dumpfile_base'] = dumpfile_base\n",
    "    results_dict['dumpfile_surf'] = dumpfile_surf\n",
    "    results_dict['E_total_base'] = E_total_base\n",
    "    results_dict['E_total_surf'] = E_total_surf\n",
    "    results_dict['A_surf'] = A_surf\n",
    "    results_dict['E_coh'] = E_total_base / system.natoms\n",
    "    results_dict['E_surf_f'] = E_surf_f\n",
    "    \n",
    "    return results_dict"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 3.3. relax_system()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "def relax_system(lammps_command, system, potential,\n",
    "                 mpi_command=None, etol=0.0, ftol=0.0, maxiter=10000,\n",
    "                 maxeval=100000, dmax=uc.set_in_units(0.01, 'angstrom'),\n",
    "                 lammps_date=None):\n",
    "    \"\"\"\n",
    "    Sets up and runs the min.in LAMMPS script for performing an energy/force\n",
    "    minimization to relax a system.\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    lammps_command :str\n",
    "        Command for running LAMMPS.\n",
    "    system : atomman.System\n",
    "        The system to perform the calculation on.\n",
    "    potential : atomman.lammps.Potential\n",
    "        The LAMMPS implemented potential to use.\n",
    "    mpi_command : str, optional\n",
    "        The MPI command for running LAMMPS in parallel.  If not given, LAMMPS\n",
    "        will run serially.\n",
    "    etol : float, optional\n",
    "        The energy tolerance for the structure minimization. This value is\n",
    "        unitless. (Default is 0.0).\n",
    "    ftol : float, optional\n",
    "        The force tolerance for the structure minimization. This value is in\n",
    "        units of force. (Default is 0.0).\n",
    "    maxiter : int, optional\n",
    "        The maximum number of minimization iterations to use (default is \n",
    "        10000).\n",
    "    maxeval : int, optional\n",
    "        The maximum number of minimization evaluations to use (default is \n",
    "        100000).\n",
    "    dmax : float, optional\n",
    "        The maximum distance in length units that any atom is allowed to relax\n",
    "        in any direction during a single minimization iteration (default is\n",
    "        0.01 Angstroms).\n",
    "    lammps_date : datetime.date or None, optional\n",
    "        The date version of the LAMMPS executable.  If None, will be identified\n",
    "        from the lammps_command (default is None).\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    dict\n",
    "        Dictionary of results consisting of keys:\n",
    "        \n",
    "        - **'logfile'** (*str*) - The name of the LAMMPS log file.\n",
    "        - **'initialdatafile'** (*str*) - The name of the LAMMPS data file\n",
    "          used to import an inital configuration.\n",
    "        - **'initialdumpfile'** (*str*) - The name of the LAMMPS dump file\n",
    "          corresponding to the inital configuration.\n",
    "        - **'finaldumpfile'** (*str*) - The name of the LAMMPS dump file\n",
    "          corresponding to the relaxed configuration.\n",
    "        - **'potentialenergy'** (*float*) - The total potential energy of\n",
    "          the relaxed system.\n",
    "    \"\"\"\n",
    "    # Build filedict if function was called from iprPy\n",
    "    try:\n",
    "        assert __name__ == pkg_name\n",
    "        calc = iprPy.load_calculation(calculation_style)\n",
    "        filedict = calc.filedict\n",
    "    except:\n",
    "        filedict = {}\n",
    "    \n",
    "    # Ensure all atoms are within the system's box\n",
    "    system.wrap()\n",
    "    \n",
    "    # Get lammps units\n",
    "    lammps_units = lmp.style.unit(potential.units)\n",
    "      \n",
    "    #Get lammps version date\n",
    "    if lammps_date is None:\n",
    "        lammps_date = lmp.checkversion(lammps_command)['date']\n",
    "    \n",
    "    # Define lammps variables\n",
    "    lammps_variables = {}\n",
    "    system_info = system.dump('atom_data', f='system.dat',\n",
    "                              potential=potential,\n",
    "                              return_pair_info=True)\n",
    "    lammps_variables['atomman_system_pair_info'] = system_info\n",
    "    lammps_variables['etol'] = etol\n",
    "    lammps_variables['ftol'] = uc.get_in_units(ftol, lammps_units['force'])\n",
    "    lammps_variables['maxiter'] = maxiter\n",
    "    lammps_variables['maxeval'] = maxeval\n",
    "    lammps_variables['dmax'] = uc.get_in_units(dmax, lammps_units['length'])\n",
    "    \n",
    "    # Set dump_modify format based on dump_modify_version\n",
    "    if lammps_date < datetime.date(2016, 8, 3):\n",
    "        lammps_variables['dump_modify_format'] = '\"%i %i %.13e %.13e %.13e %.13e\"'\n",
    "    else:\n",
    "        lammps_variables['dump_modify_format'] = 'float %.13e'\n",
    "    \n",
    "    # Write lammps input script\n",
    "    template_file = 'min.template'\n",
    "    lammps_script = 'min.in'\n",
    "    template = iprPy.tools.read_calc_file(template_file, filedict)\n",
    "    with open(lammps_script, 'w') as f:\n",
    "        f.write(iprPy.tools.filltemplate(template, lammps_variables,\n",
    "                                         '<', '>'))\n",
    "    \n",
    "    # Run LAMMPS\n",
    "    output = lmp.run(lammps_command, lammps_script, mpi_command)\n",
    "    \n",
    "    # Extract output values\n",
    "    thermo = output.simulations[-1]['thermo']\n",
    "    results = {}\n",
    "    results['logfile'] = 'log.lammps'\n",
    "    results['initialdatafile'] = 'system.dat'\n",
    "    results['initialdumpfile'] = 'atom.0'\n",
    "    results['finaldumpfile'] = 'atom.%i' % thermo.Step.values[-1]\n",
    "    results['potentialenergy'] = uc.set_in_units(thermo.PotEng.values[-1],\n",
    "                                                 lammps_units['energy'])\n",
    "    \n",
    "    return results"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4. Run calculation function(s)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "results_dict = surface_energy(lammps_command, ucell, potential, hkl,\n",
    "                              mpi_command = mpi_command,\n",
    "                              sizemults = sizemults,\n",
    "                              minwidth = minwidth,\n",
    "                              cutboxvector = cutboxvector,\n",
    "                              shiftindex = shiftindex,\n",
    "                              etol = energytolerance,\n",
    "                              ftol = forcetolerance,\n",
    "                              maxiter = maxiterations,\n",
    "                              maxeval = maxevaluations,\n",
    "                              dmax = maxatommotion)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "dict_keys(['dumpfile_base', 'dumpfile_surf', 'E_total_base', 'E_total_surf', 'A_surf', 'E_coh', 'E_surf_f'])"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "results_dict.keys()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 5. Report results"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 5.1. Define units for outputting values\n",
    "\n",
    "- __area_unit__ is the unit of area to display surface area in.\n",
    "- __energy_unit__ is the unit of energy to display cohesive energies in.\n",
    "- __energy_area_unit__ is the energy per area to report the surface energy in."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "area_unit = 'nm^2'\n",
    "energy_unit = 'eV'\n",
    "\n",
    "energy_area_unit = 'mJ/m^2'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 5.2. Print $E_{coh}$, $A_{surface}$, and $E_{surface}^f$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "E_coh =   -4.4499999983449 eV\n",
      "A_surface = 3.0975990100882553 nm^2\n",
      "E_surface_f = 1877.9483735323636 mJ/m^2\n"
     ]
    }
   ],
   "source": [
    "print('E_coh =  ', uc.get_in_units(results_dict['E_coh'], energy_unit), energy_unit)\n",
    "print('A_surface =', uc.get_in_units(results_dict['A_surf'], area_unit), area_unit)\n",
    "print('E_surface_f =', uc.get_in_units(results_dict['E_surf_f'], energy_area_unit), energy_area_unit)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
